package com.numericalmethod.algoquant.model.ralph2009.draft;

import java.io.BufferedWriter;
import java.io.FileWriter;
import java.io.IOException;
import java.text.SimpleDateFormat;
import java.util.Date;
import java.util.Random;

import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.DefaultRealMatrixChangingVisitor;
import org.apache.commons.math3.linear.RealMatrix;

public class PriceDynamics {
	protected double mu0=0;
	protected double mu1=0;
	protected double sigmaS1=0;
	protected double alpha = 0;
	protected double sigmaX1=0;
	protected double sigmaX2=0;
	protected double phi = 0; // degree of momentumness of the stock price process, range from 0 to 1
	
	protected double riskFree = 0.004;
	protected int pathCount = 100; //max: [50K, 60, 30/360]
	//protected int timePeriod = 60; //number of business days
	protected int timePeriod=3;
	protected double dt = 30.0/365; //day in year
	protected double s0 = 100; //stock price initial values, whatever values does not matter as we are interested in excess return
	protected double x0 = 0; //initial de-mean div yield assume 0 initially
	protected double m0 = 0.01; //initial momentum assume 0 initially 
	 
	
	protected double[][] z1Matrix = null; //the random number for stocks
	protected double[][] z2Matrix = null; //the 2nd components for dividend yields, since 1st component is corrected to z1
	protected RealMatrix divYieldMatrix = null;
	protected RealMatrix stockPriceMatrix = null;
	protected RealMatrix momentumMatrix = null;
	protected RealMatrix excessReturnMatrix = null;

	public static void main(String[] args){
		PriceDynamics pd = new PriceDynamics();
		pd.simulatePath();
		/*
		// !! the will throw out of memory exception, use without dumping is FINE!! 		
  		System.out.println("z1Matrix: " + pd.z1Matrix);
		System.out.println("z2Matrix: " + pd.z2Matrix);
		System.out.println("stockPriceMatrix: " + pd.stockPriceMatrix);
		System.out.println("divYieldMatrix: " + pd.divYieldMatrix);
		System.out.println("momentumMatrix: " + pd.momentumMatrix);
		*/
		
		/*
		System.out.println("z1Matrix: " + pd.z1Matrix.getColumnMatrix(pd.timePeriod-1));
		System.out.println("z2Matrix: " + pd.z2Matrix.getColumnMatrix(pd.timePeriod-1));
		System.out.println("stockPriceMatrix: " + pd.stockPriceMatrix.getColumnMatrix(pd.timePeriod-1));
		System.out.println("divYieldMatrix: " + pd.divYieldMatrix.getColumnMatrix(pd.timePeriod-1));
		System.out.println("momentumMatrix: " + pd.momentumMatrix.getColumnMatrix(pd.timePeriod-1));
		*/
	}
	
	public PriceDynamics(){
		//for testing, set to Table 1 (Panel A values for VW for Mean Reversion ONLY)
		this.phi = 0.15;
		this.mu0 = 0.92 / 100;
		this.mu1 = 0.016;
		this.sigmaS1 = 5.24 / 100;
		this.alpha = 0.011;
		this.sigmaX1 = -5.77 / 100;
		this.sigmaX2 = 1.35 / 100;
		
/*		this.phi = 0.0;
		this.mu0 = 0.15;
		this.mu1 = 0;
		this.sigmaS1 = 0.2;
		this.alpha = 0.0;
		this.sigmaX1 = 0;
		this.sigmaX2 = 0;	*/	
	}
		
	public PriceDynamics(double mu0, double mu1, double sigmaS1, double alpha,
			double sigmaX1, double sigmaX2, double phi) {
		super();
		this.mu0 = mu0;
		this.mu1 = mu1;
		this.sigmaS1 = sigmaS1;
		this.alpha = alpha;
		this.sigmaX1 = sigmaX1;
		this.sigmaX2 = sigmaX2;
		this.phi = phi;
	}
	
    public static void save(String fileName, RealMatrix matrix) {
	    try {
	
	        BufferedWriter writer = new BufferedWriter(new FileWriter(fileName));
	        for ( int i = 0; i < matrix.getRowDimension(); i++){  
	        	for(int j=0; j<matrix.getColumnDimension();j++){
	        		writer.write(matrix.getEntry(i, j) + ",");
	        	}
	        	writer.write("\n");
	        }
	        writer.close();
	    } catch(IOException ex) {
	        ex.printStackTrace();
	    }
	}
	
	public void simulatePath(){
		//assume stock price follows log normal
/*		RandomGenerator z1Rg = new JDKRandomGenerator();
		z1Rg.setSeed(17399225432l);  
		// Create a GassianRandomGenerator using rg as its source of randomness
		GaussianRandomGenerator z1RawGenerator = new GaussianRandomGenerator(z1Rg);
*/
		Random rand1=new Random(2965179875022589598L);
		
/*		RandomGenerator z2Rg = new JDKRandomGenerator();
		z2Rg.setSeed(2965179875022589598L);  
		// Create a GassianRandomGenerator using rg as its source of randomness
		GaussianRandomGenerator z2RawGenerator = new GaussianRandomGenerator(z2Rg);*/
		
		Random rand2=new Random(17399225432l);
		
		this.z1Matrix = new double[pathCount][timePeriod];
		this.z2Matrix = new double[pathCount][timePeriod];

		
		for(int i=0;i<pathCount;i++){
			for(int j=0;j<timePeriod;j++){
				z1Matrix[i][j]=rand1.nextGaussian();
				z2Matrix[i][j]=rand2.nextGaussian();
     		}
		}
		
		runDividendSimulation();
		runMomentumSimulation();
		runStockSimulation();

		SimpleDateFormat sdf=new SimpleDateFormat("yyyyMMddHHmmss");
		String dateStr=sdf.format(new Date());
		//save("C:\\dev\\MAFS\\mafs6010D\\Project\\strategic-asset-allocation\\stockPriceMatrix-"+dateStr+".csv", stockPriceMatrix);
		//save("C:\\dev\\MAFS\\mafs6010D\\Project\\strategic-asset-allocation\\divYieldMatrix-"+dateStr+".csv", divYieldMatrix);
		//save("C:\\dev\\MAFS\\mafs6010D\\Project\\strategic-asset-allocation\\momentumMatrix-"+dateStr+".csv", momentumMatrix);
		//save("C:\\dev\\MAFS\\mafs6010D\\Project\\strategic-asset-allocation\\excessReturnMatrix-"+dateStr+".csv", excessReturnMatrix);
	
	}
	
	//The S(t) simulation, i.e. Eq. 13
	// dSt = St[ ((mu0 + mu1*Xt)(1-phi) + phi*Mt)*dt + sigmaS1*dZt ] 
	protected void runStockSimulation(){
		if(stockPriceMatrix==null){
			runMomentumSimulation();
			final RealMatrix localMomentumMatrix = this.momentumMatrix;
			final RealMatrix localDivYieldMatrix = this.divYieldMatrix;
			stockPriceMatrix = new Array2DRowRealMatrix(pathCount, timePeriod);
			excessReturnMatrix = new Array2DRowRealMatrix(pathCount, timePeriod-1);
			stockPriceMatrix.walkInColumnOrder(new DefaultRealMatrixChangingVisitor(){
				public double visit(int row, int column, double value){
					if(column == 0){
						return s0;
					}else{
						double z1 = z1Matrix[row][column];
						double lastMt = localMomentumMatrix.getEntry(row, column-1);
						double lastXt = localDivYieldMatrix.getEntry(row, column-1);
						double lastSt = stockPriceMatrix.getEntry(row, column-1);
						double dSt = lastSt*( ((mu0 + mu1*lastXt)*(1-phi) + phi*lastMt)*dt + sigmaS1*Math.sqrt(dt)*z1 );
						excessReturnMatrix.setEntry(row, column-1, (dSt/lastSt) - (riskFree*dt));
						return lastSt + dSt;
					}
				}
			});
		}
	}
	
	//The M(t) simulation, i.e. Eq. 14 
	// dMt = (1 - phi)(mu0 + mu1*Xt - Mt*) *dt + sigmaS1 * dZt
	protected void runMomentumSimulation(){
		if(momentumMatrix == null){
			momentumMatrix = new Array2DRowRealMatrix(pathCount, timePeriod);
			runDividendSimulation();
			final RealMatrix localDivYieldMatrix = this.divYieldMatrix;
			momentumMatrix.walkInColumnOrder(new DefaultRealMatrixChangingVisitor(){
				public double visit(int row, int column, double value){
					if(column == 0){
						return m0;
					}else{
						double z1 = z1Matrix[row][column];
						double lastMt = momentumMatrix.getEntry(row, column-1);
						double lastXt = localDivYieldMatrix.getEntry(row, column-1);
						double dMt = (1-phi)*(mu0+ (mu1*lastXt) - lastMt) *dt + (sigmaS1 * Math.sqrt(dt) * z1);
						return lastMt + dMt;
					}
				}
			});	
		}
	}

	
	//The X(t) simulation, i.e. Eq. 15 (Xt follows OU process)
	// dXt = -alpha * Xt * dt + sigmaX' * dZt and 
	// mut = mu0 + mu1 * Xt
	// or, 
	// Xt+1 = Xt + dXt
	protected void runDividendSimulation(){
		if(divYieldMatrix == null){
			divYieldMatrix = new Array2DRowRealMatrix(pathCount, timePeriod);
			divYieldMatrix.walkInColumnOrder(new DefaultRealMatrixChangingVisitor(){
				public double visit(int row, int column, double value){
					if(column == 0){
						return x0;
					}else{
						double z1 = z1Matrix[row][column];
						double z2 = z2Matrix[row][column];
						double lastXt = divYieldMatrix.getEntry(row, column-1);
						double dXt = (-alpha * dt * lastXt) + (sigmaX1*Math.sqrt(dt)*z1) +(sigmaX2*Math.sqrt(dt)*z2);
						return lastXt + dXt;
					}
				}
			});		
		}
	}

}
